home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / int_3d.pro < prev    next >
Text File  |  1997-07-08  |  12KB  |  258 lines
  1. ;$Id: int_3d.pro,v 1.7 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       INT_3D
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes the triple integral of a trivariate
  11. ;       function F(x,y,z) with limits of integration from A to B
  12. ;       for X, from P(x) to Q(x) for Y and from U(x,y) to V(x,y)
  13. ;       for Z.
  14. ;
  15. ; CATEGORY:
  16. ;       Numerical Analysis.
  17. ;
  18. ; CALLING SEQUENCE:
  19. ;       Result = INT_3D(Fxyz, AB_Limits, PQ_Limits, UV_Limits, Pts)
  20. ;
  21. ; INPUTS:
  22. ;       Fxyz:  A scalar string specifying the name of a user-supplied
  23. ;              IDL function that defines the trivariate function to be
  24. ;              integrated. The function must accept x, y & z and return
  25. ;              a scalar result.
  26. ;       AB_Limits:  A two-element vector containing the lower, A, and
  27. ;                   upper, B,  limits of integration for x.
  28. ;       PQ_Limits:  A scalar string specifying the name of a user-
  29. ;                   supplied IDL function that defines the lower, P(x),
  30. ;                   and upper, Q(x), limits of integration for y. The
  31. ;                   function must accept x and return a two-element
  32. ;                   vector result.
  33. ;       UV_Limits:  A scalar string specifying the name of a user-
  34. ;                   supplied IDL function that defines the lower, U(x,y),
  35. ;                   and upper, V(x,y), limits of integration for z. The
  36. ;                   function must accept x & y and return a two-element
  37. ;                   vector result.
  38. ;        Pts:  The number of transformation points used in the
  39. ;              computation. Possible values are: 6, 10, 20, 48 or 96.
  40. ;
  41. ; KEYWORD PARAMETERS:
  42. ;       DOUBLE: If set to a non-zero value, computations are done in
  43. ;               double precision arithmetic.
  44. ;
  45. ; EXAMPLE:
  46. ;       Compute the triple integral of the trivariate function
  47. ;       F(x,y,z) = z * (x^2 + y^2 + z^2)^1.5 over the region:
  48. ;       A = -2, B = 2, Px = -sqrt(4 - x^2), Qx = sqrt(4 - x^2),
  49. ;       Uxy = 0, Vxy = sqrt(4 - x^2 - y^2).
  50. ;      
  51. ;       ;Define the trivariate function.
  52. ;         function Fxyz, x, y, z
  53. ;           return, z * (x^2 + y^2 + z^2)^1.5
  54. ;         end
  55. ;
  56. ;       ;Define the limits of integration for y.
  57. ;         function PQ_Limits, x
  58. ;           return, [-sqrt(4. - x^2), sqrt(4. - x^2)]
  59. ;         end
  60. ;
  61. ;       ;Define the limits of integration for z.
  62. ;         function UV_Limits, x, y
  63. ;           return, [0.0, sqrt(4. - x^2 - y^2)]
  64. ;         end
  65. ;
  66. ;       ;Define the limits of integration for x.
  67. ;         AB_Limits = [-2.0, 2.0]
  68. ;
  69. ;       ;Integrate with 10, 20, 48, and 96 point formulas using double-
  70. ;       ;precision arithmetic. Notice that it is possible to abbreviate
  71. ;       ;keywords.
  72. ;         print, INT_3D('Fxyz', AB_Limits, 'PQ_Limits', 'UV_Limits', 10, /d)
  73. ;         print, INT_3D('Fxyz', AB_Limits, 'PQ_Limits', 'UV_Limits', 20, /d)
  74. ;         print, INT_3D('Fxyz', AB_Limits, 'PQ_Limits', 'UV_Limits', 48, /d)
  75. ;         print, INT_3D('Fxyz', AB_Limits, 'PQ_Limits', 'UV_Limits', 96, /d)
  76. ;
  77. ;       INT_3D with 10 transformation points yields:    57.444248
  78. ;       INT_3D with 20 transformation points yields:    57.446201
  79. ;       INT_3D with 48 transformation points yields:    57.446265
  80. ;       INT_3D with 96 transformation points yields:    57.446266
  81. ;       The exact solution (6 decimal accuracy) yields: 57.446267
  82. ;
  83. ; PROCEDURE:
  84. ;       INT_3D.PRO computes the triple integral of a trivariate function
  85. ;       using iterated Gaussian Quadrature. The algorithm's transformation
  86. ;       data is provided in tabulated form with 15 decimal accuracy.
  87. ;
  88. ; REFERENCE:
  89. ;       Handbook of Mathematical Functions
  90. ;       U.S. Department of Commerce
  91. ;       Applied Mathematics Series 55
  92. ;
  93. ; MODIFICATION HISTORY:
  94. ;       Written by:  GGS, RSI, January 1994
  95. ;       Modified:    GGS, RSI, September 1994
  96. ;                    Added 96 point transformation data.
  97. ;                    Added DOUBLE keyword. Replaced nested FOR loop with
  98. ;                    vector operations resulting in faster execution.
  99. ;       Modified:    GGS, RSI, April 1996
  100. ;                    Modified keyword checking and use of double precision.
  101. ;-
  102.  
  103. FUNCTION Int_3D, Fxyz, AB_Limits, PQ_Limits, UV_Limits, Pts, Double = Double
  104.  
  105.   ON_ERROR, 2
  106.  
  107.   if N_ELEMENTS(AB_Limits) ne 2 then $
  108.     MESSAGE, "AB_Limits parameter must be a two-element vector."
  109.  
  110. ; Tabulated transformation data with 15 decimal accuracy.
  111. if Pts eq 6 then begin
  112.   Ri    = DBLARR(Pts)          &   Wi    = DBLARR(Pts)
  113.   Ri[0] = 0.932469514203152d   &   Wi[0] = 0.171324492379170d
  114.   Ri[1] = 0.661209386466265d   &   Wi[1] = 0.360761573048139d
  115.   Ri[2] = 0.238619186083197d   &   Wi[2] = 0.467913934572691d
  116.   Ri[INDGEN(Pts/2) + (Pts/2)] = - Ri[(Pts/2) - INDGEN(Pts/2) -1]
  117.   Wi[INDGEN(Pts/2) + (Pts/2)] =   Wi[(Pts/2) - INDGEN(Pts/2) -1]
  118. endif else if Pts eq 10 then begin
  119.   Ri    = DBLARR(Pts)          &   Wi    = DBLARR(Pts)
  120.   Ri[0] = 0.973906528517172d   &   Wi[0] = 0.066671344308688d
  121.   Ri[1] = 0.865063366688985d   &   Wi[1] = 0.149451349150581d
  122.   Ri[2] = 0.679409568299024d   &   Wi[2] = 0.219086362515982d
  123.   Ri[3] = 0.433395394129247d   &   Wi[3] = 0.269266719309996d
  124.   Ri[4] = 0.148874338981631d   &   Wi[4] = 0.295524224714753d
  125.   Ri[INDGEN(Pts/2) + (Pts/2)] = - Ri[(Pts/2) - INDGEN(Pts/2) -1]
  126.   Wi[INDGEN(Pts/2) + (Pts/2)] =   Wi[(Pts/2) - INDGEN(Pts/2) -1]
  127. endif else if Pts eq 20 then begin
  128.   Ri     = DBLARR(Pts)         &   Wi     = DBLARR(Pts)
  129.   Ri[0]  = 0.993128599185094d  &   Wi[0]  = 0.017614007139152d
  130.   Ri[1]  = 0.963971927277913d  &   Wi[1]  = 0.040601429800386d
  131.   Ri[2]  = 0.912234428251325d  &   Wi[2]  = 0.062672048334109d
  132.   Ri[3]  = 0.839116971822218d  &   Wi[3]  = 0.083276741576704d
  133.   Ri[4]  = 0.746331906460150d  &   Wi[4]  = 0.101930119817240d
  134.   Ri[5]  = 0.636053680726515d  &   Wi[5]  = 0.118194531961518d
  135.   Ri[6]  = 0.510867001950827d  &   Wi[6]  = 0.131688638449176d
  136.   Ri[7]  = 0.373706088715419d  &   Wi[7]  = 0.142096109318382d
  137.   Ri[8]  = 0.227785851141645d  &   Wi[8]  = 0.149172986472603d
  138.   Ri[9]  = 0.076526521133497d  &   Wi[9]  = 0.152753387130725d
  139.   Ri[INDGEN(Pts/2) + (Pts/2)] = - Ri[(Pts/2) - INDGEN(Pts/2) -1]
  140.   Wi[INDGEN(Pts/2) + (Pts/2)] =   Wi[(Pts/2) - INDGEN(Pts/2) -1]
  141. endif else if Pts eq 48 then begin
  142.   Ri     = DBLARR(Pts)         &   Wi     = DBLARR(Pts)
  143.   Ri[0]  = 0.998771007252426d  &   Wi[0]  = 0.003153346052305d
  144.   Ri[1]  = 0.993530172266350d  &   Wi[1]  = 0.007327553901276d
  145.   Ri[2]  = 0.984124583722826d  &   Wi[2]  = 0.011477234579234d
  146.   Ri[3]  = 0.970591592546247d  &   Wi[3]  = 0.015579315722943d
  147.   Ri[4]  = 0.952987703160430d  &   Wi[4]  = 0.019616160457355d
  148.   Ri[5]  = 0.931386690706554d  &   Wi[5]  = 0.023570760839324d
  149.   Ri[6]  = 0.905879136715569d  &   Wi[6]  = 0.027426509708356d
  150.   Ri[7]  = 0.876572020274247d  &   Wi[7]  = 0.031167227832798d
  151.   Ri[8]  = 0.843588261624393d  &   Wi[8]  = 0.034777222564770d
  152.   Ri[9]  = 0.807066204029442d  &   Wi[9]  = 0.038241351065830d
  153.   Ri[10] = 0.767159032515740d  &   Wi[10] = 0.041545082943464d
  154.   Ri[11] = 0.724034130923814d  &   Wi[11] = 0.044674560856694d
  155.   Ri[12] = 0.677872379632663d  &   Wi[12] = 0.047616658492490d
  156.   Ri[13] = 0.628867396776513d  &   Wi[13] = 0.050359035553854d
  157.   Ri[14] = 0.577224726083972d  &   Wi[14] = 0.052890189485193d
  158.   Ri[15] = 0.523160974722233d  &   Wi[15] = 0.055199503699984d
  159.   Ri[16] = 0.466902904750958d  &   Wi[16] = 0.057277292100403d
  160.   Ri[17] = 0.408686481990716d  &   Wi[17] = 0.059114839698395d
  161.   Ri[18] = 0.348755886292160d  &   Wi[18] = 0.060704439165893d
  162.   Ri[19] = 0.287362487355455d  &   Wi[19] = 0.062039423159892d
  163.   Ri[20] = 0.224763790394689d  &   Wi[20] = 0.063114192286254d
  164.   Ri[21] = 0.161222356068891d  &   Wi[21] = 0.063924238584648d
  165.   Ri[22] = 0.097004699209462d  &   Wi[22] = 0.064466164435950d
  166.   Ri[23] = 0.032380170962869d  &   Wi[23] = 0.064737696812683d
  167.   Ri[INDGEN(Pts/2) + (Pts/2)] = - Ri[(Pts/2) - INDGEN(Pts/2) -1]
  168.   Wi[INDGEN(Pts/2) + (Pts/2)] =   Wi[(Pts/2) - INDGEN(Pts/2) -1]
  169. endif else if Pts eq 96 then begin
  170.   Ri     = DBLARR(Pts)         &   Wi     = DBLARR(Pts)
  171.   Ri[0]  = 0.999689503883230d  &   Wi[0]  = 0.000796792065552d
  172.   Ri[1]  = 0.998364375863181d  &   Wi[1]  = 0.001853960788946d
  173.   Ri[2]  = 0.995981842987209d  &   Wi[2]  = 0.002910731817934d
  174.   Ri[3]  = 0.992543900323762d  &   Wi[3]  = 0.003964554338444d
  175.   Ri[4]  = 0.988054126329623d  &   Wi[4]  = 0.005014202742927d
  176.   Ri[5]  = 0.982517263563014d  &   Wi[5]  = 0.006058545504235d
  177.   Ri[6]  = 0.975939174585136d  &   Wi[6]  = 0.007096470791153d
  178.   Ri[7]  = 0.968326828463264d  &   Wi[7]  = 0.008126876925698d
  179.   Ri[8]  = 0.959688291448742d  &   Wi[8]  = 0.009148671230783d
  180.   Ri[9]  = 0.950032717784437d  &   Wi[9]  = 0.010160770535008d
  181.   Ri[10] = 0.939370339752755d  &   Wi[10] = 0.011162102099838d
  182.   Ri[11] = 0.927712456722308d  &   Wi[11] = 0.012151604671088d
  183.   Ri[12] = 0.915071423120898d  &   Wi[12] = 0.013128229566961d
  184.   Ri[13] = 0.901460635315852d  &   Wi[13] = 0.014090941772314d
  185.   Ri[14] = 0.886894517402420d  &   Wi[14] = 0.015038721026994d
  186.   Ri[15] = 0.871388505909296d  &   Wi[15] = 0.015970562902562d
  187.   Ri[16] = 0.854959033434601d  &   Wi[16] = 0.016885479864245d
  188.   Ri[17] = 0.837623511228187d  &   Wi[17] = 0.017782502316045d
  189.   Ri[18] = 0.819400310737931d  &   Wi[18] = 0.018660679627411d
  190.   Ri[19] = 0.800308744139140d  &   Wi[19] = 0.019519081140145d
  191.   Ri[20] = 0.780369043867433d  &   Wi[20] = 0.020356797154333d
  192.   Ri[21] = 0.759602341176647d  &   Wi[21] = 0.021172939892191d
  193.   Ri[22] = 0.738030643744400d  &   Wi[22] = 0.021966644438744d
  194.   Ri[23] = 0.715676812348967d  &   Wi[23] = 0.022737069658329d
  195.   Ri[24] = 0.692564536642171d  &   Wi[24] = 0.023483399085926d
  196.   Ri[25] = 0.668718310043916d  &   Wi[25] = 0.024204841792364d
  197.   Ri[26] = 0.644163403784967d  &   Wi[26] = 0.024900633222483d
  198.   Ri[27] = 0.618925840125468d  &   Wi[27] = 0.025570036005349d
  199.   Ri[28] = 0.593032364777572d  &   Wi[28] = 0.026212340735672d
  200.   Ri[29] = 0.566510418561397d  &   Wi[29] = 0.026826866725591d
  201.   Ri[30] = 0.539388108324357d  &   Wi[30] = 0.027412962726029d
  202.   Ri[31] = 0.511694177154667d  &   Wi[31] = 0.027970007616848d
  203.   Ri[32] = 0.483457973920596d  &   Wi[32] = 0.028497411065085d
  204.   Ri[33] = 0.454709422167743d  &   Wi[33] = 0.028994614150555d
  205.   Ri[34] = 0.425478988407300d  &   Wi[34] = 0.029461089958167d
  206.   Ri[35] = 0.395797649828908d  &   Wi[35] = 0.029896344136328d
  207.   Ri[36] = 0.365696861472313d  &   Wi[36] = 0.030299915420827d
  208.   Ri[37] = 0.335208522892625d  &   Wi[37] = 0.030671376123669d
  209.   Ri[38] = 0.304364944354496d  &   Wi[38] = 0.031010332586313d
  210.   Ri[39] = 0.273198812591049d  &   Wi[39] = 0.031316425596861d
  211.   Ri[40] = 0.241743156163840d  &   Wi[40] = 0.031589330770727d
  212.   Ri[41] = 0.210031310460567d  &   Wi[41] = 0.031828758894411d
  213.   Ri[42] = 0.178096882367618d  &   Wi[42] = 0.032034456231992d
  214.   Ri[43] = 0.145973714654896d  &   Wi[43] = 0.032206204794030d
  215.   Ri[44] = 0.113695850110665d  &   Wi[44] = 0.032343822568575d
  216.   Ri[45] = 0.081297495464425d  &   Wi[45] = 0.032447163714064d
  217.   Ri[46] = 0.048812985136049d  &   Wi[46] = 0.032516118713868d
  218.   Ri[47] = 0.016276744849602d  &   Wi[47] = 0.032550614492363d
  219.   Ri[INDGEN(Pts/2) + (Pts/2)] = - Ri[(Pts/2) - INDGEN(Pts/2) -1]
  220.   Wi[INDGEN(Pts/2) + (Pts/2)] =   Wi[(Pts/2) - INDGEN(Pts/2) -1]
  221. endif else MESSAGE, "Pts parameter must be 6, 10, 20, 48 or 96."
  222.  
  223.   TypeAB = SIZE(AB_Limits)
  224.  
  225.   ;If the DOUBLE keyword is not set then the internal precision and
  226.   ;result are identical to the type of input.
  227.   if N_ELEMENTS(Double) eq 0 then $
  228.     Double = (TypeAB[TypeAB[0]+1] eq 5) 
  229.  
  230.   if Double eq 0 then begin
  231.     Ri = FLOAT(Ri) & Wi = FLOAT(Wi)
  232.   endif
  233.  
  234.   H1 = (AB_Limits[1] - AB_Limits[0])/2.0
  235.   H2 = (AB_Limits[1] + AB_Limits[0])/2.0
  236.   Aj = 0.0
  237.   for i = 0, Pts-1 do begin
  238.     X  = H1 * Ri[i] + H2
  239.     jX = 0.0
  240.     CF = CALL_FUNCTION(PQ_Limits, X)
  241.     k1 = (CF[1] - CF[0])/2.0
  242.     k2 = (CF[1] + CF[0])/2.0
  243.     for j = 0, Pts-1 do begin
  244.       Y  = k1 * Ri[j] + k2
  245.       CF = CALL_FUNCTION(UV_Limits, X, Y)
  246.       l1 = (CF[1] - CF[0]) / 2.0
  247.       l2 = (CF[1] + CF[0]) / 2.0
  248.       jY = TOTAL(Wi * CALL_FUNCTION(Fxyz, X, Y, l1*Ri+l2), Double = Double)
  249.       jX = jX + Wi[j] * l1 * jY
  250.     endfor
  251.     Aj = Aj + Wi[i] * k1 * jX
  252.   endfor
  253.  
  254.   if Double eq 0 then RETURN, FLOAT(Aj * H1) else $
  255.     RETURN, (Aj * H1)
  256.  
  257. end
  258.